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I. INTRODUCTION 

In recent years there has been a growth in the under¬ 
standing of the importance of intrinsic noise in com¬ 
plex systems composed of a finite number of interact¬ 
ing constituents. It has been recognised that there 
are visible, macroscopic effects due to the stochastic- 
ity inherent in the interactions in a variety of sys¬ 
tems, including for example metabolic pathways [1], 
gene regulatory systems ms, predator-prey dynam¬ 
ics HIT] or models of disease spread [SI IS] ■ This in¬ 
herent stochasticity is referred to as ‘intrinsic noise’ 
or ‘demographic stochasticity’ m- Intrinsic noise 
can give rise to phenomena such as cyclic dynamics 
[3 E], patterns and waves HMS], and extinction 
events dlHlI]- These phenomena are not captured by 
more traditional deterministic modelling approaches, 
instead they are purely noise-induced. 

Deterministic models are built on ordinary or par¬ 
tial differential equations. These equations are for¬ 
mally only valid in the limit of an infinite system, 
that is the number of interacting constituents is so 
large that stochastic effects play no role [22] • To 
take into acount the intrinsic stochasticity of the in¬ 
teractions a full probabilistic description is required. 
Widely used modelling approaches are drawn from 
the theory of stochastic processes, most notably the 
master equation [22l [23] , describing the time evolu¬ 
tion of the probability distribution over the space of 
states. Formulating the master equation approach 
relies on the Markov property of the underlying dy¬ 
namics: transition rates from one state to another 
must only depend on the target state and the present 
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state of the system, but not on the path the dynamics 
has taken to arrive at the current state. This implies 
that the system has no memory of previous interac¬ 
tions, all effects of an interaction must be realised 
instantaneously. While this is a reasonable assump¬ 
tion for many processes in physics, this is typically 
not the case in biological models. For example birth 
in a predator-prey system occurs after a period of 
gestation, transcriptional and translational delays are 
relevant in gene regulatory systems |24j . and recov¬ 
ery occurs a certain period of time after infection in 
models of disease spread [25l|26]. These processes can 
lead to situations in which effects of an initial inter¬ 
action (e.g. impregnation, initiation of transcription, 
infection) materialise with a significant delay. Such 
dynamics are then no longer Markovian as the change 
of the state of the system at any one time t may de¬ 
pend on what processes were set in motion in the 
past and which complete at t. Of course the mod¬ 
elling as a delay system can often be avoided through 
the construction of higher-dimensional models with 
intermediate states (e.g. staged models in epidemics 
[26]), but mathematically it is often convenient to use 
delay models, as these are relatively easy to set up for 
arbitrary delay distributions. 

Traditionally delay dynamics have been investigated 
based on deterministic approaches [MlElHSIl. These 
are subject to the limitations outlined above in that 
they do not capture the effects of stochasticity. It 
is only recently that analytical (and also numeri¬ 
cal) techniques for individual-based models subject 
to both intrinsic noise and delay have been devel¬ 
oped . In previous work we have derived Gaus¬ 

sian approximations of such dynamics, and we have 
shown that these can capture the effects of delay re¬ 
actions gni SI]. This analysis is conveniently car¬ 
ried out in terms of generating functionals [131 [E] , 
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and the equivalent of the linear-noise approximation 
(LNA) for delay systems can be derived. An alterna¬ 
tive, more informal method to derive the same Gaus¬ 
sian approximations (by-passing the exact descrip¬ 
tion) was highlighted in [iT] . 

In these existing approaches delay reactions fire at 
an initial time with rates determined by the state 
of the system at that time. The reaction may then 
have an immediate effect on the composition of the 
system, and subsequently a second effect materialises 
at a later time (after the delay). The delay is either 
fixed, or it can be drawn from an underlying distribu¬ 
tion each time a delay reaction is initiated. Existing 
work is often restricted to what we will refer to as 
‘definitive’ completion of the delay reaction. Once a 
delay reaction is initiated (it ‘fires’) the delayed ef¬ 
fect will occur, no matter what the trajectory of the 
system is between the time the reaction fires and the 
designated completion time. This is an obvious lim¬ 
itation of the modelling approach. Descriptions in 
which delay reactions may fail to complete depend¬ 
ing on events that occur between initiation and des¬ 
ignated completion provide a much more realistic de¬ 
scription of many real-world processes. The gestation 
period in a predator-prey model presents perhaps the 
most intuitive example of delay reactions which can 
be interrupted. A pregnancy is initiated at a given 
time and the birth event will occur at a designated 
later time. However, birth is not certain, the mother 
might die during pregnancy. Similarly, in a model of 
disease spread, an individual may get infected with 
the disease and would then be scheduled to recover at 
a later time. However, the individual may die in the 
interim, or be removed from the system in some other 
way, so that the completion event (recovery) may not 
occur. The rate with which such removal occurs may 
well depend on the state of the system at the time of 
removal. 

The purpose of the present paper is to extend existing 
approaches to the modelling of stochastic dynamics 
with delay to the case in which delay reactions may 
not complete. We will refer to such reactions as ‘in¬ 
terruptible delay reactions’. In this setting delay re¬ 
actions can be of several types: (i) Delay reactions 
which cannot be interrupted. This is the case consid¬ 
ered in 00]; ( ii) Delay reactions which can be inter¬ 
rupted, but the probability of interruption does not 
depend on the state of the system. An example is 
death of the mother for internal reasons. Another 
example was also considered in m in the context of 
the susceptible-infective-recovered (SIR) model with 
birth and death, with constant death rate; (iii) Delay 
reactions which can be interrupted and the probabil¬ 
ity of interruption depends on the state of the system 
at the time of interruption. Continuing the example 
of a pregnant prey individual, the rate of predation 
depends on the number of predators present in the 
system throughout the pregnancy period. This case 


is not covered by [10] . 

In this paper we develop a systematic Gaussian ap¬ 
proximation to models with interruptible delay re¬ 
actions. Specifically we extend the generating func¬ 
tional method used in |40j to include interruption ef¬ 
fects so that all three cases above are covered. As 
an application we consider the SIR model of disease 
spread and two variants of a predator-prey model, one 
with a delay period due to pregnancy and one with a 
delay period due to the maturation of juveniles. 


II. MODEL DEFINITIONS: 

INTERRUPTIBLE DELAY REACTIONS 

II.1. Model definition 

We consider a finite population of interacting con¬ 
stituents, each of which is of one of M different types, 
labelled a = We will refer to the con¬ 

stituents as ‘individuals’ in the following. The state 
of the system at any one time is then described by 
an M-dimensional vector, n(t) = [ni{t),... ,nM(t)], 
where the non-negative integer na(t) indicates the 
number of individuals of type a at time t. We as¬ 
sume a continuous-time evolution. The system is 
well-mixed and two individuals of the same type are 
indistinguishable. The individuals interact through a 
set of R reactions, we will label these i = 1,... ,R. 
Each possible reaction i is associated with an initia¬ 
tion rate T((n), indicating the rate with which reac¬ 
tions of type i fire if the system is in state n. When 
such a reaction fires an instantaneous change of the 
state of the system n{t) occurs at the time of firing, 
described by the vector v^. That is to say the state of 
the system changes from n to n -|- . Subsequently 

these reactions may also have a delayed effect. This 
is implemented as follows: if a delay reaction fires at 
time t an instantaneous change of the state of the sys¬ 
tem occurs as described above. In addition to this, a 
delay time t is drawn from an underlying delay dis¬ 
tribution, Ki{T). This distribution is specific to the 
reaction triggered, as indicated by the subscript i. 
After the delay time has elapsed a second change in 
the state of the population may occur. This happens 
at time t + t, and we denote the delayed effect of 
the reaction by . The key element of the dynamics 
that we add in the present work is the possibility that 
a delay reaction may not complete. We assume that 
a delay reaction of type i, triggered at time t and due 
to complete at t + t can be interrupted at any time 
between t and t + r. The rate with which this occurs 
is fi[n{t')], where t <t' < t + r. The termination rate 
may thus depend on the state of the system n(<'). If 
this happens the delayed effect at time t + r does not 
occur, instead we assume that the state of the system 
changes by at time t'. For the time being we will 
only consider cases in which there is only one way 
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in which each delay reaction can be interrupted. A 
generalisation to models in which delay reactions can 
be interrupted in multiple different ways is relatively 
straightforward, though slightly more cumbersome, 
see Appendix [B] for details. 

As is commonly done in the modelling of interact¬ 
ing particle systems, we will assume that all reaction 
rates Ti (n) scale with a parameter - that is to say 
Ti(n) = 0(fl). This parameter can be seen as setting 
the scale of the size of the system (the scale of the 
total number of individuals in the system, or equiv¬ 
alently the volume of the system), and the scaling 
of the rates reflects the fact that the total (average) 
number of reactions in the system per unit time is 
proportional to its size. For later purposes it is con¬ 
venient to introduce the quantities Xa{t) = na(t)/fl. 
We will refer to these as the ‘concentrations’ of parti¬ 
cles of type a. We also introduce the intensive rates 
r^[x{t)] = Ti[x{t)]/n. 


time f (0 < s < t). We note that non-delay reactions 
are included in this descriptions, they would simply 
have Wi = u,; = 0. 

Eq. (§ defines the conditional probability 


P{xt+A\{xt,k,£, m}t'<t) 

= n <5 Ut+A.a - ^ ^ ^ 


i r>0 




0<s<r 


(3) 


where the notation {x^k,£,m}ti<t indicates vari¬ 
ables with indices t' <t. 

We have already described the Poissonian nature of 
the variables ki^t = tiut it remains to define in 

more detail how the and are chosen. 

This is explained in the following section. 


II.2. Discrete-time dynamics 


II.3. Statistics of reaction numbers 


As in [40j we will proceed by discretising time into 
steps of duration A. The continuum limit will be 
restored at the end. We will write x^^t for the con¬ 
centration of individuals of type a at time step t. In 
the discretised model, we will assume that all reac¬ 
tion rates remain constant between t and < -I- A. If 
the concentration vector is Xt at time t then the num¬ 
ber of newly triggered reactions of type i during this 
time step is a Poissonian random variable ki^t with 
parameter Ti(n)A. In the absence of delay reactions 
the dynamics of the discrete-time model would hence 
read 




( 1 ) 


In models with delay reactions we have to extend 
this expression to take into account (i) delay reac¬ 
tions completing at the designated time and (ii) de¬ 
lay reactions which terminate before the designated 
completion time. We then have 


^i+A,a ^t,a. — q EE 


T>0 


Xt,aklt + 


We first discuss the statistics of the variables kj^, in¬ 
dicating the number of reactions of type i triggered in 
the discrete time model at time step t and with com¬ 
pletion due at time t + t. As discussed in [ID] each 
is a Poissonian random variable with parameter 
A^Ki{T)£lri{xt). This is not affected by possible in¬ 
terruptions of the delay reactions. At this point it 
is useful to recall the definition ri{xt) = Ti{xt)/£l. 
Broadly speaking £lri{xt) reflects the rate with which 
reactions of type i are initiated (with any delay). 
Once a reaction is initiated a delay period is drawn 
independently from a distribution Kilr). In the 
discrete-time model this is reflected in the factor 
AKi{T). 

Ultimately we will be interested in taking the 
continuous-time limit for the dynamics. Anticipating 
this, we focus on the case of small A, which simplifies 
the problem significantly. The probability distribu¬ 
tion of kl^ is of the form 

P{kl, = l\xt) = A^K,{T)nn{xt) + 0{A'^), 

P{klt = 0|a;t) = 1 - A^K,{T)nr,{xt) + ©(A^), 
Pikl, > l\x^) = 0{A^). (4) 


+ Y 

0<S<T 


( 2 ) 


In this expression kj^ is the number of reactions of 
type i that fire at time t and which have a delay 
period t. The quantity is the number of reac¬ 

tions of type i that fired at time t — r and successfully 
complete at time t. Finally, the quantity is the 
number of reaction of type i that hred at time t — s 
with a designated delay period r (i.e. scheduled for 
completion at t — s + r), but which are interrupted at 


Eq. Q indicates that the probability to observe two 
or more initiation events of reactions of type i and 
with delay r in any one time interval A is at least of 
order A"^. In the limit of small A we can therefore 
restrict ourselves to kj^ € {0,1}. 

If kl^ = 0 then it is clear that = 0 for all s and 
also = 0, i.e. we have 

P(£[;/=0|fc:,, = 0) = lV0<s<r, 

P{ml, = Q\klt = 0) = 1. (5) 
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If = 1 then only one out of the {^^’/}o<s<r and 
can be non-zero. If kl^ = I and = 1 then it 
follows that = 0 for all s ^ a, and m\^ = 0. I.e., 
we have 


their statistical weight. Knowing the weight of each 
combination enables us to average over the fc, and 
m; (iv) Finally, the continuous-time limit is restored. 
The resulting generating functional is 


m:: = m:: = i, = i) = i v a < s < r, 

= = l,kl = l) = l. (6) 


= J DxDp ( 11 ) 


As explained above, interruption happens with prob¬ 
ability A X fi{xt+s) in the next time interval A if the 
system is in state Xt+s, and so 

P{£l’J = l\xt+,, = 0}o<.<s, kl^ = 1) 

P{£i: = 0\xt+,, {tl’J = 0}o<.<«, kl, = 1) 

= l-/\f,{xt+,). (7) 

Finally, if the reaction has not been interrupted by 
time t + T then the reaction always completes, 

P{ml, = = 0}o<.<., kl = 1) = 1. (8) 


III. GENERATING FUNGTIONAL 

In discrete time the generating function is 

^(-0) = , (9) 

\ / paths 

where the average is performed over all possible 
paths. We have introduced a source term 0, deriva¬ 
tives with respect to it generate correlation functions 
[48l |44] The main task in the calculation that follows 
is to compute the path-average in the above expres¬ 
sion. This entails integrating/summing over all ran¬ 
dom variables, 

Z(0) = ^ 

t kX,m 

( 10 ) 

where V{x, k,i,m) is the joint probability distribu¬ 
tion of all Xt^a, , and - i.e. the probability 

of a path. Our objective is to find an expression for 
the generating functional after performing the sums 
over all k, £, and m and after subsequently taking 
the continuous-time limit A —>■ 0. The algebra in¬ 
volved is somewhat lengthy, and so we do not give 
full details here, they are relegated to Appendix fA| 
The calculation consists of carrying out the following 
main steps: (i) The path probability V{x, k,£,m) is 
expressed as a product of conditional probabilities us¬ 
ing Eqs. (ii) The delta-functions in Eq. @ 

are converted into their Fourier representations, in 
the process introducing the conjugate variables p*; 
(iii) Using Eqs. the combinations of variables 

with non-zero probability are identified, along with 


with the action 


^[a;, p] = —y^t 




( 2 ) 


( 12 ) 


There are two contributions to the action from each 
type of reaction: (i) a contribution from the reactions 
which are not interrupted, 

- Rl^\t) = y 

X e- ko d^M^(t-<^)^Ki{T)nr,[x{t - t)], (13) 

and (ii) a contribution from the reactions which are 
interrupted, 

- i?f)(t) = y dsJ^e-AK p(*-«)+“. p(*)l _ 

pOO 

xMx{t)]e-ko‘i-M-(t-.)] / dTK,^T)nn[x{t - s)]. 

J s 

(14) 


These contributions are both functionals of the path 
of X between the time at which reaction first fires and 
the time it either completes successfully or is inter¬ 
rupted. 

We notice the factor e~ ko d<Tfi[x{t-a)] (13). 


This exponential is the probability that a delay re¬ 
action of type i triggered at t — r and with com¬ 
pletion time t, reaches completion, given a path 
X between t — t and t. Similarly, the quantity 
f,[x{t)]e- -fo [cf. Eq. indicates the 

probability that a delay reaction of type i triggered at 
time t — s is interrupted in the time interval [t, -|-fit), 
given a path x between t — s and t. 

If all delay reactions of type i complete with certainty 
(i.e. in absence of interruptions) one has fi{x) = 0. 
This implies =0 (see Eq. (El)- The quantity 


(t) in Eq. (13) reduces to the corresponding ob¬ 
ject in 


IV. APPLICATION TO THE SIR MODEL 
WITH BIRTH AND DEATH 

The SIR model with birth and death can be studied 
using the approach developed in Sec. |n] This model 
describes the dynamics of an infectious disease in a 
population of individuals. Each individual can be in 
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one of three states: susceptible, infectious or recov¬ 
ered. We consider an infection dynamics with delayed 
recovery. Upon contact with an infectious individual 
a susceptible member of the population may become 
infectious (with rate /3), and then they recover at a 
time T after infection. The delay time t is drawn 
from a recovery-time distribution K(t). We write 
this process as follows 

S + I ^21-1^ R. (15) 

The first arrow indicates the effect the reaction has 
at the time it is initiated. The double arrow indicates 
the delayed effect r units of time after the reaction 
is triggered. In addition to the recovery process, any 
individual may die, this occurs at constant rate /x, 
independent of the infection status of the individual. 
In order to keep the population size, N, constant any 
dying individual is immediately replaced by an indi¬ 
vidual of type S. This leads to the following addi¬ 
tional reactions 

I ^S, 

S. (16) 

Crucially, an individual may die (and be replaced by 
an S) after becoming infected, but before the sched¬ 
uled recovery time. This represents a delay reaction 
with possible interruption in the formalism described 
in the previous sections. In this introductory exam¬ 
ple, however, the rate of interruption does not de¬ 
pend on the state of the system and so this model 
can be studied using the simpler method presented 
m- There we mapped the above reaction scheme 
onto the model 


R-!^S, 

S + I R, 


with 


(17) 


/*oo poo 

x = J dTK{T) J ds 

/ oo 

ds /xe"'"®, 

pOO 

Q(s) = (1 - x)"Ve"^® / dr K{t). 


(18) 


The first line in Eq. 0 describes the standard 
death and replacement reaction R S. The sec¬ 
ond and third lines represent two reaction sequences 
which may be triggered when an infection event oc¬ 
curs. Each sequence starts with an infection event 
(occurring with rate /3njns/N). With probability x 
the newly infected individual is scheduled for recovery 
at a later time; the time-to-recovery, r, is drawn from 


the distribution K{t). With complementary proba¬ 
bility 1 — X the newly infected individual is scheduled 
for death (and replacement by an individual of type 
S). The time-to-replacement, s, is drawn from the 
distribution Q{s). This mapping is possible because 
the interruption rate, /x, is independent of the state 
of the system. Both distributions, K and Q, are nor¬ 
malised. 


The reaction scheme in Eq. 0 and the distribu¬ 
tions in Eq. (181 were formulated in [30] , and from 
this scheme the generating functional of the process 
was derived. Now that we have put the more gen¬ 
eral formalism of the previous sections in place we 
can recognise these existing results as a special case. 
Inserting the model specifications into the general for¬ 
malism of the previous sections the action in Eq. (12) 
is indeed found as 


S[x, P] = y dti^Ps(t)is(t) + pi{t)xi{t) 

+ — l^fV(l — xs{t) — xi{t)) 

+ J dtye"wh(‘')-ps(t')-px(t)] _ 

xK{t - t')Nx/3xs{t')xi(t') 

+ J dt'w[pid')-psit')-piit)+ps{t)] _ 

xQ{t - t')N{l - x)Pxs{t')xi{t') i. 


(19) 


where we have used the definitions of Xi and 

Q{s) as given in Eqs. (18). The total size of the pop¬ 
ulation, N = ns -\- ni + nn is constant in time, so 
there are only two independent degrees of freedom. 
In formulating Eq. (191 we have therefore elimi¬ 
nated species R via nji = N — ns — nj (equivalently 
xr = 1 — xs — xi). It is important to stress that we 
have used U = iV for simplicity, given the constant 
size of the population. The result in Eq. (19), derived 
from the general formalism of the previous sections, 
is the same action as the one found in |3D] using the 
mapping of Eq. 0. 


V. APPROXIMATIONS TO THE 
GENERATING FUNCTIONAL 

V.l. Deterministic approximation 


The exact generating functional can be approximated 
by Taylor expanding the action in powers of the in¬ 
verse system size. To leading order in 12“^ the general 
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action Eq. (12) reduces to 


S'[a:, p] = — J dt 
where 




( 20 ) 


v4f ^ (t) = j dr • p{t -t) +Wi- p{t) 

xe-fJ - r)] 


( 21 ) 


and 


Af ^ (t) = J ds • p{t - s) +Ui- p{t) 

x/i[a;(t)]e--/'o 

poo 

X / dTKi{T)ri[x{t - s)]. (22) 

J S 

This action is linear in p (by construction), and it 
describes the deterministic dynamics 

= Y.v.,^n[x^{t)]+FW[x^]{t)+Fi^^[x^]it), 

(23) 


with the drift terms 


_ r poo 

EW[a;“](t)=V / dr 

. Jo 


X K,{T)ri[x°°{t - r)] 


, (24) 


and 


r poo 

Fi^)[x^]{t) = Y^ / dsuxcf[x^{t)] 

, Ldo 

pOO 

X g-/o d^/dx”(t-a)] / dTK,{T)n[x°°{t - S)] 

J s 

(25) 


featuring due to the delay. The deterministic approx¬ 
imation is valid for large (formally infinite) system 
sizes, and neglects all stochastic effects. We have de¬ 
noted the dynamical variables in the deterministic 
limit by x°°. 


V.2. Linear-noise approximation 


variables via p —>■ The resulting generating 

functional is 


Z[i/)] = J D^Dq 


(26) 


Differentiating Eq. (26) with respect to ipit) still gives 
the moments and correlation functions of x(t). To 
find a generating functional for ^ we divide Eq. ( |26[ ) 
by the factor e~ ^-ip(t)-x (t)^ which is independent 
of ^ and q, and rescale the source term, ^l} ^/Ucj). 

The generating functional for £ is 


Z^[<p] = [ D^Dq ^^^C.vT2q]-/di</)(t)-$(t)^ 

. 

The moments of ^ can now be found by differentiating 
Z^[cj)\ with respect to 4>- 

The terms in the action are again expanded in powers 
of and the expansion is curtailed after sub¬ 

leading order (i.e. keeping terms 0(12°) and above). 
In order to illustrate the procedure we will consider 
the case where the deterministic dynamics are at a 
stable fixed point, x°°{t) = x*. This is not nec¬ 
essary for the linear-noise approximation, however 
it will allow us to keep the amount of algebra un¬ 
der control. The action is S'[a:* -1- -\/f2q] = 

-S'LNAil, q] + 0(f2“l/2) 


= - f dt'^qa{t)ia{t) 

Ol 



. In this expres- 


with Ac,^p[x*\{t,t') = 
sion stands for the right-hand side of Eq. (|23|). 


The action is now quadratic in q and and it de¬ 
scribes a process with additive Gaussian coloured 
noise [HI 05] . This process is of the form 



dt' Aa^l3[x*]{tA’)ip{t') +qa{t), 


(29) 


The linear-noise approximation is used to study fluc¬ 
tuations about the solution to the deterministic equa¬ 
tions of motion, Eq. ( |2^ . We write x = a;°°-|-f2“^/^^, 
and formulate the generating functional in terms of 
the variable In doing this we re-scale the conjugate 


with noise correlator {qa{t)qp[t')) = Ba^p[x*]{t,t')- 
The explicit expressions for Aa^d[x*]{t,t') and 
Ba,p[x*\{tA') are somewhat lengthy, and before pre¬ 
senting them it is helpful to make a few simplifying 
definitions. 













7 


V.3. Further simplications 


To simplify the above expressions further we define 
the distributions 


Qi(s) = 


1 * 

dTK,(T), 


(30) 


for T,s > 0. We set Ki(T) = 0 for t < 0, 
and Qi(s) = 0 for s < 0, and we introduce 
Xi{x*) = dTe~^d^ )TX-(^r). Within the 
linear-noise approximation and assuming that the 
deterministic dynamics are at the fixed point x* 
this is the probability that a delay reaction of type 
i, once triggered, completes. The function Ki{T) is 
the conditional time to completion. The quantity 
1 — Xi{x*), on the other hand is the probability that 
the delay reaction is interrupted before completion, 
and Qiir) describes the conditional distribution 
of times to interruption. Both distributions Ki{T) 
and Qiir) are normalised. One notices that the 


expressions in Eq. (18) are special cases of Eq. (30) 


with fi{x*) = /r. We stress though that our analysis 
below will explicitly account for (linear) fluctuations 
of the interruption rate as the state of the system 
varies in time. These fluctuations of the interruption 
rate are obviously not present in the system with 
constant fi{x) = considered in [iO] . 


Additionally it is helpful to make the definition 

Li,air) = Xi{'^*)Ki{T)Wi,a 

+[l - Xi{x*)]Qi{T)ui,o,. (31) 


This quantity can be interpreted as the weighted de¬ 
lay effects (due to either successful completion or in¬ 
terruption) on species a of reactions of type i a time 
period r after such a reaction was triggered. 

Using these definitions we can write Aa,p[x*\{t,t’) 
concisely as 




+Ui,a[l - Xi{x*)] 

POO 

- / dr Li,air) 


dri{x*) 

S{t t) 
ri{x*) dfi{x*) 


Mx*) dx*^ 
dft{x*) 


S{t -1') 


dxt 




r^{x*)Q{t - t') 


(32) 


Each of the four terms in the curly brackets has a 
clear interpretation. As the system fluctuates around 
the fixed point, the number of reactions triggered, in¬ 
terrupted and completed will fluctuate as well. The 


first term is the contribution to ^a(t) of initial ef¬ 
fects of reactions (at the time of triggering) due to 
fluctuations. The second term is the contribution of 
interruption effects due to fluctuations at the time 
of interruption. The third term is the contribution 
of delayed effects (either successful completion or in¬ 
terruption) due to fluctuations between the reaction 
firing and the delayed effects occurring. Finally, the 
fourth term is the contribution of delayed effects due 
to fluctuations at the time the reaction first fired. 


If the rate of interruption, fi, is independent of the 
state of the system for all i then only the first and 
last terms in Eq. (32) are non-zero. 


We now turn to the correlation Ba,i 3 [x*](t,t') of the 
Gaussian noise variables in the linear-noise approx¬ 
imation. Given that delay reactions can have ef¬ 
fects on the composition of the population at multiple 
times, this noise is not white. Instead one finds 


Ba,i3[x*]{t,t') = ^ri(x*)<^ V^,aVi,p + Wi,aWi,ljXi{x*) 

i ^ 

+U^,aU^,0[l - Xlix*)] 

+Vi,aLi,p{t' -t)+ Vi,pLi,a{t - t')|- (33) 


As Li,a{t) = 0 for t < 0 only one of the three terms 
in Eq. (33) is non-zero for any pair of times t and L 
and any given reaction i. 


V.4. Spectrum of fluctuations 


The linear-noise approximation can be used to calcu¬ 
late the Fourier spectrum of fluctuations about the 
deterministic fixed point [7] . This is particularly use¬ 
ful to study noise-induced quasi cycles, as we will 
explain in further detail below. Fourier transforming 
Eq. (29), gives 


iw^Q(w) = X! 1 y,aXi(x*) - Li,a(uj) 

+ ( g.o + L,,^{uj)) 

+ ^ |?;3(^) + ^a(w). 

(34) 


This equation is linear in and it can be written as 
fj{uj) = M(w)4(a;), with a suitable matrix M(a;). The 
Fourier transform of Eq. (33) (with respect to t — t') 



















is 

Ba,p{^) = V^^aVi^p + Wi,aWi,/3Xi(®*) 

i ^ 

“t“ W2,a^z,/3[1 X*(^ )] 

(35) 

The spectrum of fluctuations about the deterministic 
fixed point is then characterised by the matrix S(w) = 

and it can be found from (see [46]) 
S(w) =M(a;)-i]B(w)[Ml'(w)]-^ (36) 


VI. APPLICATION TO A PREDATOR-PREY 
MODEL 


VI.1. Model definitions 

As an application of the above formalism we will now 
consider stochastic effects in a predator-prey model 
with different types of delay. Specifically we write 
X to denote prey-individuals, and Y for predators. 
We focus on the dynamics governed by the following 
reactions 


X 2A, 

Y + X 2Y, 

Y 0. (37) 

We write nx and ny for the number of individuals of 
each type, and x = nxji^,y = ny/ft. As before fl is 
a parameter, setting the scale of the population size. 
The first reaction describes reproduction of prey with 
a logistic birth rate, dependent on the concentration 
X. The constant k represents a carrying capacity, 
more precisely, the system can contain a most fcO in¬ 
dividuals of the prey type. The logistic birth rate 
for prey distinguishes this model from other stochas¬ 
tic predator-prey models [7]. Within our stylised ap¬ 
proach we assume that predation events result in the 
birth of a predator (second reaction), the third reac¬ 
tion finally describes a death process for predators. 
This is obviously a minimalist model, but it is in line 
with previous stylised modelling approaches for birth- 
death processes ST], and it serves as a helpful testbed. 
There are various ways in which delay processes can 
be introduced into this model. One is to include 
gestation periods, in which a prey enters a pregnant 
state before giving birth. After the gestation period 
the pregnant individual returns to the regular (non¬ 
pregnant) state and a new prey individual is created. 
Another possibility is to assume that newly born prey 
individuals are initially in a ‘juvenile’ state and that 


they cannot immediately reproduce. After a matu¬ 
ration period they become full adult prey individuals 
and acquire the ability to reproduce. We will refer 
to the first modification as the ‘gestation model’ and 
the second modification as the ‘juvenile model’. In 
both models there is an additional intermediary class 
of individuals, denoted X' - pregnant prey in the ges¬ 
tation model and juvenile prey in the juvenile model. 
The class X corresponds to non-pregnant prey in the 
gestation model and adult prey in the juvenile model. 
We write nx' for the number of individuals of type 
X' and and x' = nx> /fl. Of course one could intro¬ 
duce a similar state in the context of predators, but 
in-line with the above stylised approach we keep the 
complexity of the model to a minimum. 

For the gestation model the reactions are 


x'.x'^2X, 

Y + X ^2Y, 

Y -^0. 


(38) 


where we have introduced Xtot = x + x', the total 
concentration of prey. In the first reaction, a non¬ 
pregnant individual becomes pregnant with rate 6(1 — 
Xtot/k). They then give birth after a delay drawn 
from K(t) to a non-pregnant individual and return 
to the non-pregnant state. 

For the juvenile model we use 




Y + X -^2Y, 
Y ^0. 


(39) 


In this model an adult prey gives birth to a juvenile 
individual with rate 6(1—cctot /k) . After a delay drawn 
from K{t) the juvenile becomes an adult individual. 
Additionally, in both models, the intermediary indi¬ 
viduals, X', can be predated upon, via 

Y + X' ^ 2Y. (40) 


It is this predation reaction that leads to the possibil¬ 
ity of interrupting a delay reaction before its sched¬ 
uled completion. If a pregnant individual is elimi¬ 
nated in a predation event, it will obviously not give 
birth at the designated time, and similarly if a juve¬ 
nile individual is removed during predation it will not 
reach its reproductive age. Crucially, the rate with 
which such events happen depends on the state of the 
system (e.g. on the number of predators in the pop¬ 
ulation at that time). The methods of [10] are hence 
not applicable, and we will instead use the extended 
and more general formalism developed earlier in this 
paper. 

Both models have the same reaction rates. The only 
difference between the two models lies in the changes 
of nx and nx' when a delay reaction triggers and 
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completes, i.e. in the numerical values of the stoi¬ 
chiometric coefficients vi^x and wi^x (we label the de¬ 
lay reaction by i = 1 in both models). It is therefore 
convenient to carry out the analysis for both mod¬ 
els together, keeping vi^x and wi^x general. After 
performing the calculations we will substitute in for 
these parameters and compare the two models. 
Simulations of the models can be performed using 
the Modified Next Reaction Method (MNRM) 
[34l [35] . The interruption reaction shown in Eq. (401 
fires like a conventional reaction in the MNRM 
algorithm with rate fi[x{t)] x mi{t) where mi{t) 
is the number of ‘active’ delay reactions at time t, 
i.e. reactions of type i = 1 which have fired but 
for which the delayed effects have not yet occurred. 
When the interruption reaction occurs an element 
of the fist of queued delay reactions is selected at 
random with uniform weights, and removed from 
the list. The effects of the interruption reaction 
are then applied to the state of the system accord¬ 
ing to ui^a- Aside from this the MNRM is unchanged. 


VI.2. Deterministic dynamics 


Both models describe three species, A, A' and Y. 
The primary reactions are listed in Eq. (38) and (391 
respectively. These are labelled i = 1, 2, 3 from top 
to bottom. The reactions rates for both systems are 


/ x{t)h[xtot{t)] \ 

r[x{t)] = px{t)y{t) , (41) 

V Mt) j 


with h{xtot) = b{l — Xtot/k). The first reaction, i = 
1, is the only delay reaction. This reaction can be 
interrupted due to predation on the intermediaries, 
with rate fi[x{t)] = py(t) per intermediary (i.e. any 
instance of the delay reaction which is active at time t 
is subject to interruption with rate py{t)). The delay 
distribution is K{t). The stoichiometric coefficients, 
Vi^ai describing changes to the system when reactions 
trigger can be summarised as follows. 


( Vx I 0 \ 

-10 1 , (42) 

VO 0 -ly/ 


where the rows each stand for one reaction (i = 
1,2,3) and the columns represent the three types of 
individuals (A, A' and Y). The only non-zero sto- 
chiometric coefficients for interruption events [see Eq. 
(40)] are ui^xtat — Ui^y = 1. The non-zero 

stochiometric coefficients for successful completion of 
delay reactions are wi^x = Wx and wi^x' = ~1- The 
differences between the two models are in the numer¬ 
ical values of Vi^x = W and Wi^x = Wx- Eor the juve¬ 
nile model we have Vx = O-iWx = 1, for the gestation 
model we have Vx = —1, Wx = 2. 


The deterministic equations of motion are found as 


a;°°(t) = VxX°°{t)h[x'^^{t)]-px°°{t)y°°{t) 

pOC) 

+Wx / dre- ^0 d<ypy°°it-a) 

Jo 

xK{T)x^{t-T)h[x^^{t-T)], (43) 


y^{t)=px°^{t)y°^{t)-dy^{t) 

nOO 

+ / dr py°°{t)e-^o 

Jo 

/ oo 

dsK{s)x°°{t - T)h[x^^^{t - t)], (44) 

with = x°° + x'°°. In principle we can also write 
an equation for x '^, however we know from the reac¬ 
tions that the number of X' at any one time t is equal 
to the total number of A' created up to that point 
(through birth) and which have not been removed 
through maturation or predation. We find 


X 


/OO 



e~ fo d<ypy°°(t-(7) 


ds K{s)x^{t-T)h[x^^{t 



VI.3. Fixed point analysis 


At the fixed point the left-hand sides of Eqs. (43) 
and (44) are zero, the concentration of intermediary 
individuals can be found from Eq. (li^. We obtain 


0 =[dJxx{y*) + - px*y*, (46) 

0 = px*toty* - dy*, (47) 

0 = [1 - x(2/*)]a:*/i(a;tot) - Px'*y*, (48) 


where the stars indicate the fixed-point values of the 
relevant variables. We have also introduced the quan¬ 
tity x{y*) = dr e~Py’'^K{T), representing the 
probability for a delay reaction to reach completion 
when the system is at the fixed point. Eor the gesta¬ 
tion model this is the fraction of pregnant individuals 
which successfully give birth, for the juvenile model 
this is the fraction of juveniles which successfully ma¬ 
ture and reach the reproductive age. 

There are three solutions to the above fixed point 
equations: all extinct {x = Xtot = 2 / = a:' = 0), only 
prey (x = Xtot = k, y = x' = 0), and coexistence. 
We will focus on the coexistence fixed point, and we 
will simply refer to it as the fixed point. We see 
that = d/p for both models and for any delay 
distribution. The coexistence fixed point only exists 
if d/p < k. If d/p > k then hid/p) = (1 — ^) < 0 
and Eqs. (46)-(|4^ have no consistent solutions. If 
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d/p < k we can rearrange Eqs. (46)-(48), and find 

(49) 

(50) 


y*=p ^[w^xiv* 


■Vx]h{-), 


X* = (j) [wxxiy*) + Vx]-, 

x'* = (j)-^[l-x{y*)\-, 

P 


(51) 


with (j) = vjxxi.y*) + Vx + 1 — x{y*)- We will now 
restrict the further discussion to the case of constant 
delay f, i.e. K{t) = — f). In this case we can 

proceed with the analysis, and find x{y*) = 

Eqs. (49)-([M|) are a transcendental set of equations 
for X *, x'* and y* , which cannot easily be simpli¬ 
fied any further. However we can reparametrise the 
model, and treat x as a parameter in Eqs. (49)-(51). 
The delay period, r, is then a function of the model 
parameters through f = —hi{x)/(jpy*)- After substi¬ 
tuting in for y* we find 


r = - ln(x)[(w„x + Vx)Kp\ 


d\^-l 


(52) 



FIG. 1. Relationship between survival probability at the 
fixed point, x, and delay period, r. Lines are parametric 
plots using Eq. (52l. Blue lines and circles correspond 


to the gestation model, yellow lines and triangles to the 
juvenile model. Parameters are: b = p = k = 1 and 
d = 0.2. For the gestation model the survival proba¬ 
bility approaches 0.5, whereas for the juvenile model it 
approaches 0. Symbols are the fraction of delay reactions 
which completed in simulations using the MNRM, aver¬ 
aged over 100 realisations and with 12 = 1000. 


If X = 1 then there is no delay (f = 0) in either model. 
Furthermore x decreases with increasing delay r. In 
the gestation model x 0.5 as t —>• oo, whereas 
X —t 0 in the juvenile model, as shown in Fig. We 
stress that we have not formally established stability 
of fixed point. However for the parameters used in 
Fig.[T] the fixed point has numerically been seen to 
be stable up to at least f = 5. 

In Fig. we compare these theoretical predictions 
against data from simulations of the microscopic dy¬ 
namics in finite populations (12 = 1000). As seen in 
the figure the above deterministic analysis is in good 
agreement with stochastic simulations of the micro¬ 
scopic process in a finite system, 12 = 1000. 

The concentration of predator and prey individuals 
X* and y*, determined by Eqs. (49)-([^ decrease 


with the delay period in both models, whereas the 
concentration of intermediary individuals, x'* , natu¬ 
rally increases when the delay period becomes longer. 
This is summarised in Fig. obtained as a paramet¬ 
ric plot of Eqs. (|4^-([^ and Eq. ([52|. Comparison 


with exact stochastic simulations at finite 12 again 
shows good agreement. We notice that the concen¬ 
tration of predators is much more sensitive to the 
choice of model than the concentration of prey for 
non-zero delay. 


VI.4. Linear-noise approximation 

The above deterministic analysis is valid only in the 
limit of infinite populations. We next study the ef¬ 
fects of noise induced by the stochastic dynamics 
when the number of individuals in the system is finite. 
Representative trajectories generated from stochastic 
simulations in Fig.j^show that the effects of noise can 



FIG. 2. Relationship between fixed point location and 
delay peri od, r. Li nes are parametric plots using Eq. ( |52| | 
and Eqs. (49l-(51l. Blue lines and circles correspond to 
the gestation model, yellow lines and triangles to the juve¬ 
nile model. Parameters are: b = p = k — 1 and d = 0.2. 
The total prey concentration, is not affectected by 
the delay (black line), and is the same for both models. 
Symbols are from simulations using the MNRM averaged 
over 100 realisations with XI = 1000. 


be quite profound. As seen in numerous Markovian 
systems [HlllIillllllH] noise can generate sustained 
oscillations in parameter regimes in which the purely 
deterministic system approaches a stable fixed point. 
This effect has also been observed in stochastic dy¬ 
namics with fixed and distributed delay [311 133MT] , 
the time series shown in Fig. [^show that this phe¬ 
nomenon extends to delay dynamics with uncertain 
completion of delay reactions. 

The oscillations can be characterised by the power 
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spectra of fluctuations about the deterministic fixed 
point, more precisely by the diagonal elements od 
S(a;) = introduced before Eq. (36). 

Fig. a shows how the power spectrum responds to 
an increase in the delay period. In both models the 
effect of delay is to increase the period of the quasi¬ 
cycles, as can be seen by a shift of the peak towards 
lower frequencies with increasing delay. Delay also 
increases the amplitude of the oscillations. For the 
gestation model the effect is more pronounced, the 
amplitude is much more sensitive to the delay. 

To calculate the power spectrum of fluctuations an¬ 
alytically we start from the general expression for 
the Fourier transform of the Langevin equation found 
in the linear-noise approximation, Eq. (34). For the 
predator-prey models we have 



-Py*) + Tx(w))^x(w) 

bx* ( ~ N — 

- px* (^1 -f Ga,{uj)ly{uj) + Tjxiuj), 

(53) 

iw'ioo'iuj) = {h{xlot) - (l + Lx' 

-+ Lj;l (W)^ ^x' (t) 

- (px* + px*h{x*^^)^Gx'{u})^y{t) + rj^'{u}), 

(54) 

—' / bx^ N ~ ^ bx^ ~ —' 

iuj^yiuj) = [py* - —jLy{uj)^^[uj) - —Ly{uj)^x'iuj) 

- px*h{xl^^)Gy{uj)^y{t) + 7jy{u}), (55) 


where La{uj) = dr e “'^^^(r) and Ga{oj) = 
^ /o°° dr Lair) {l — For these models the 

weighted delay effects are L{t) = [wxK{t),—{K{t)-\- 
Q{t)),Q{t)\. 


Eqs. (53)-(^55b can be rewritten as a matrix equation 


riijjj) = M(a;)4(a;). The noise correlation matrix has 
elements 


Bx,x{uj) =c(vl+ xwl + WxX 

+Vx + 2vxWxReK , 

Bx,x'i(^) =c(vx- XWx +w^K{uj) 

-VxiK*{uj) + Q*{uj))y 

Bx,y{iXl) = C (vxQ*iuj) - WxX - Vx'j , 

B,',x'{w) = 2c (f - Re(i?(w) + Q(w))) , 
Bx',y{uj) = C (Q*{uj) - (1 - X)) , 

By^y{uj) = 2c(j>, (56) 


FIG. 3. Sample trajectories of x using the gestation model 
for T = 0 (top line) and r = 5 (bottom line). Parameters 
are: 6 = p = fc = 1 and d = 0.2. Simulations performed 
using the MNRM with = 10000. Both trajectories are 
observed to oscillate about their respective deterministic 
fixed points. The oscillations for f = 5 are much more 
defined, with a period T ~ 50. 


where c = x*h[xl^^). The remaining elements fol¬ 


low from the Hermitian property, ]B(a;) = BT(a;). 
The spectrum matrix S(a;) can be found by inserting 
Eqs. ([^-([^ into Eq. (1^. Comparison of these 


theoretical predictions against data from simulation 
shows good agreement, see Fig. and confirms the 
viablity of the linear-noise approximation for both 
models. 

If we return back to the definition of the models, the 
delay reaction for the gestation model is 




b{l -xt^t /k) ^ K{t ) 


2X, 


whereas for the juvenile model it is 






X. 


(57) 


(58) 


The total effect of the delay reaction (i.e. of the initial 
effects together with the delayed effects) is the same 
for both models: the number of individuals of type X 
increases by one. The difference between the models 
is in the intermediary changes. 

To probe the effects of delay, let us compare the out¬ 
come of the delay models (t > 0) with the corre¬ 
sponding Markovian models (obtained in the limit 
f —> 0). In this limit the reactions in Eq. (57) and 
Eq. (1^ both reduce to 




b(^ l — x^ k) 


2X, 


(59) 


The study of the fixed point behaviour and of the 
power spectrum of the models with delay shows that 
the location of the fixed point and the peak of the 
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CO 


FIG. 4. Variation of the power spectrum of x, ^(a;) with 
the delay period, r. Top panel corresponds to the gesta¬ 
tion model, bottom panel to the juvenile model. Dashed 
lines are theoretical predictions using Eq. (521, solid lines 
are from simulations using the MNRM with = 10000. 
Parameters are: b — p = k = 1 and d = 0.2. 


power spectrum at any fixed f > 0 are closer to 
the corresponding quantities at f = 0 in the juvenile 
model than in the gestation model. This indicates 
that it is not only the duration of the delay which 
determines whether or not it can be neglected, but 
also the details of the delay reaction. We note that 
the survival probability at the fixed point (shown in 
Fig. 0 is more sensitive to the delay in the juvenile 
model than in the gestation model. 


VII. CONCLUSIONS 

The main focus of this paper has been to extend the 
generating-functional approach for stochastic inter¬ 
acting particle systems to dynamics in which delay re¬ 
actions can be interrupted. In particular we consider 
cases in which the interruption rate depends on the 


state of the system at the time of interruption. The 
probability with which the delay effects ultimately 
occur is then a functional of the path taken by the 
system during the delay period. We successfully set 
up a path-integral description of such systems. 

The action of the generating functional for delay reac¬ 
tions with interruption has a clear relationship with 
the actions found previously in models in which de¬ 
lay reaction complete with certainty. The example 
in Sec. |IV| shows how our extended theory includes 
results previously derived in simpler cases where the 
interruption rate is constant and independent of the 
state of the system at the time of interruption. The 
generating functional provides a starting point for 
further analytical calculations, specifically we derive 
Gaussian and linear-noise approximations from it to 
characterise noise-induced effects in finite popula¬ 
tions. 

We apply these techniques to an example inspired 
by predator-prey dynamics. We study two variants 
of the model, one with gestation delay and one with 
maturation delay. The deterministic approximation 
can be used to give a good estimate of the probabil¬ 
ity that any given delay reaction completes. In the 
gestation model this describes the probability that a 
pregnant individual successfully gives birth, in the ju¬ 
venile model it is the probability that a newborn indi¬ 
vidual successfully reaches adulthood. Starting from 
the linear-noise approximation we show how the for¬ 
malism can be used to analytically study persistent 
noisy cycles, induced by demographic stochasticity. 
Finally, the overall effect of delay reactions is iden¬ 
tical in the gestation and in the juvenile predator- 
prey model, that is to say the changes of delay re¬ 
actions at the time they trigger taken together with 
the effects at completion are the same in both mod¬ 
els. Our analysis shows that different quantities (e.g. 
survival probability, fixed point location, peak of the 
power spectrum) have different sensitivities to the de¬ 
lay in the two models. For both the deterministic 
fixed point and the characteristic frequency of noise- 
induced oscillations, the duration of the delay period 
has a much stronger effect on the outcome in the ges¬ 
tation model than in the juvenile model. Conversely, 
the probability that the delayed effects successfully 
occur decreases much faster with increasing delay in 
the juvenile model than in the gestation model. 
These observations have implications for the con¬ 
struction of population models. Whether a delay can 
safely be neglected depends not just on the duration 
of the delay period itself but also on the details of the 
changes to the population at the start and end of the 
delay period. This indicates that care must be taken 
when approximating models with delay or intermedi¬ 
ate processes by effective Markovian dynamics. 

In summary, our previous work |40| together with the 
present paper provides a systematic and coherent ap¬ 
proach to studying stochastic effects in a wide range 
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of delay systems, including those with distributed de¬ 
lay, delay reactions that can fail to complete with 
path-dependent rates and systems with multiple pos¬ 
sible interruption channels. The generating func¬ 
tional is the natural mathematical entity which with 
to address these systems, and is the analog of the 
master equation in the context of Markovian dynam¬ 
ics. We have here only applied this method to a small 
set of examples, inspired by ecological dynamics, but 
we anticipate that the formalism that is now in place 


can be of value in the context of many other appli¬ 
cations, including those in epidemiology, metabolic 
dynamics, gene expression and other fields which de¬ 
layed dynamics. 
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Appendix A: The generating functional for interruptible delay 


1. Conditional probabilities 


To derive the generating functional we focus on the model in discretised time. The continuous-time limit is 
taken at the end of the calculation. We adopt the following notation: kt denotes all kl^, for which t' = t, £t 
denotes all for which f + s = t, and rrit denotes all for which t' + t = t. 

Using this notation the equation of motion for Xt^a can be expressed in the form Xt+A,a = 9a{kt,£t, ^^t)- 

The path probability P{x, k,£, m) can be written as the product 

P{x, k, £, m) = n /c, -6, 7Tl\ , mf\Xt , {x, /c, TTl} t'<t) • (^1) 

t 


For a particular all of the variables contained in m^} are independent of each other (conditioned on 

the history of the system up to time t), so their joint probability distribution factorises 


, TTi/: I , {x, A;, TH-<^) — nn 


i T>0 ^ 0<S<T 

All of these conditional probabilities are known from the definition of the process. We have 

P{mlt_^\Xt,{x,k,£,m}t'<t) = 

P{£l;t-s\^t, {x, fc, £, = P{£l’t_Jxt, {£l;t_,}o<a<s,klt_,) 

P{kL\xt,{x,k,£,m}t,<t) = P{kL\xt). 


t'<t) 


(A2) 


(A3) 


The explicit expressions for the probabilities in terms of fi{xt), ri{xt) and ^^(r) are given by Eqs. 

The probability distribution P(a;(_|_A|{a:, fc, £, is the product of ^-distributions given in Eq. which 

for clarity we repeat here, 


P{xt+A\{x,k,£,m}t'<t) = P[dfx4+A,a - - 9aikt,£t,mt)y 

rv N / 


We have 


{kt,£t,mt) = Vi^akJ^t + + Y 


n 




i r>0 0<s<r 

To make the expressions that follow more compact we can define the product 

= p{kiAxt)Pimy\{ei’;}o<.^r,ky) n Pie:’;\xt+,,{eino<.<s,kit), 


0<S<T 


SO that the path probability factorises into 


P{x, fc, £, m) = \lP{xt+A\{x,k,£,m}t'<t) n 

t T>0,i,t 


(A4) 


(A5) 


(A6) 


(A7) 
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We will refer to ^ht as the statistical weight associated with the combination of random variables 
}o<cr<T) It is not itself a probability, as the probabilities which compose it are conditioned on 

different random variables. 


All of the conditional probabilities which make up in Eq. (A6| are only conditioned on variables which 
are either i) {xt’}t<t><t+T or ii) other members of {kJ^, The weight is therefore only a 

function of and {xti}t<t><t+T- To simplify the notation we suppress these arguments. 

Later on this dehnition of will allow us to factorise the generating functional and perform the sums for 

each factor separately. 


Inspecting Eqs. ©-(HI) we can identify the only combinations of values of kl^, 
is non-zero. They are as follows: 


|T,0<S<T 
t ’ 


and mj f for which 


(i) A reaction of type i with delay period r does not fire at time t. In this case we have 


Kt = 


= 0V0<o-<r, 


mlt = 0> 


with 


= 1 - A'^Ki{T)nr^{xt). 


(A8) 

(A9) 


(ii) A reaction of type i with delay period r hres at time t and is interrupted at time t + s. The values of the 
random variables for this case are 

= 0V0<CT<s 

C;; = 1, 

= 0Vs<(T<r, 

"lit = 0. 


and 


t r = 


Afi{Xt+s)llo<a<s (l - Afi{xt+a)'^A^Ki{T)nn{xt). 


(AlO) 

(All) 


(iii) A reaction of type i with delay period r fires at time t and is not interrupted. The random variables take 
the values 




= 0V0<o-<r, 


"il.t = 


with 




{l-Afi{xt+a)^A^Ki{T)Vtri{xt). 


0<<7<T 


(A12) 

(A13) 


2. Generating functional 


Inserting Eq. (A7) and Eq. (|A4|) into the definition of the generating functional gives 


i r>0 *- 


= Dx E n EE + E 

{k,£,m} t,Oi ^ 

n ‘f’i.t.T 


0<s<r 




r>0,2,t 


(AM) 
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where we have introduced the short-hand J Dx = Using the exponential representation of the 

(^-function we can write 


/ VxDx 6xp j 

*E**.“i 


^ t,Oi ^ 


^i+A,a ^t,o 






i r>0 


0<s<r 


t,ct 


(A15) 


with/f?* = /nti^. 

to give 


In what follows we will make the transformation ix — >■ p. Eq. (A15l can be factorised 


E n ^ 

{k,e,m} T>0,i,t 


-AEo 




xe 


-nEo 


Eo<«<^“..-Pt+<. 




(A16) 


Before we take the average over the random variables k, I, and m it is useful to define 


'kM.. = EE E 

UU {U’f }o<s<x 


-nE. 


Vi,aPt,akJ^+Wi^aPt,am1^ 


-nEo 


Eo<s<x “i,cPt + s.aU,t 


(A17) 


The quantity is a function of {xt'}t<<t'<t+T, for clarity we have suppressed these arguments. In defining 

we have collected all occurrences of the variables indicated in the summation, {kli,{il’f}o<cr<T,m'lf}. 
As explained beneath Eq. (A7), is only a function of {£p('^}o<cr<r, wj’d and of {a:t'}t<<t'<t+T. This 
construction allows us to evaluate each separately to all other terms in the generating functional, which 

can be written as 


Z[iP] = JVxDp 


(A18) 




Using Eqs. (A9), (All), and (A13) we obtain 


Cx = 1 


E 

0<S<T 


o- h Ec, K,cPt,o-|-Mi,cPt+s,c] 


A/,(a;t+^) {I -/^fi{xt+a)) 


0<cr<s 


0<(7<T ^ 


(A19) 


which for small A can also be written as an exponential, cf. 1 -I- A/ = + 0{A‘^). 

Repeating this procedure for all leaves a generating functional for the discrete-time dynamics. We can 

then take the continuous-time limit, and arrive at 




(A20) 


with the action 


S[X, p] = Jdt p{t) ■ X{t)+ Y[J^ g- 


+ ds p(‘)I _ fi{x{t))e- daU{a^(t-a)) ^ dr Ki{T)VLr,{x{t - s)) 


• (A21) 
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Appendix B: Multiple interruption reactions 


Although the notation is involves a number of indices, it is straightforward to see how the generating functional 
can be extended to the case in which any delay reaction can be interrupted in multiple different ways. If reaction 
f is a delay reaction with which can be interrupted in possible ways, indexed by /r = 1,..., A^, then by 
extension of Eq. ([^ , 


Xt+A 




i r>0 


Vi,aklt + 




(Bl) 


We also have to modify the conditional probabilities of the random variables to take into account the different 
interruption reactions. If kl^ = 0 then ^ = 0 for all s and /x, i.e. 


= 0) = 1 V 0 < s < r , /X = 1,..., A,. 

Similarly, if kj^ = 1 and = 1 then necessarily = 0 and = 0, i.e. 

= 0|C:;,t = 1, fcl,, = 1) = 1 V a < s < r , XX e A„ 

^Kx = oic:;,t = i:fcM = i) = i- 

If the system is in state Xt then interruption through channel /x happens with rate 

= 0}o«x<s. u=l,...MAlt = 1) = A/i,^(Xt+s), 

Pill’ll,t ~ Ol^^t+s, ~ 0}o<cr<s, ixGAi, = I) = 1 — A/i_^(a:t_|_s). 

If the reaction has not been interrupted then the reaction always completes, 

PilTT-lt = = 0}o«T<r, u=l,...,Ai,kJ^t = I) = 1. 

As we are working with small A, 

pie?Y = 1. {C,A.t = o}o«.<x, x=i,...,A.,ki, = 1) = <5^... 

The only non-zero combinations of these conditional probabilities are: 

(i) Reaction does not fire: 

Pimlt = m,t = 0) n = ^)PiKt = 0|a;t) = 1 - A^K,iT)Qn{xt). 

0<s<r 

(ii) Reaction fires and is interrupted at time f -f s through channel /x: 

Pimlt = kl = 1) n = 1’ Kt = 1) 

s<a<T 

v=l,...,Ai 

^ Pi^Ali,t ~ l|®t+s) = 0}o<o-<s, i/=l,...,Ai) fciy = I) 

= 0}o<p<.. = ^)Piklt = l|a^x) 


0<cr<s 

v£Ai 


Afi^^{xt+s) n (l - A/x.,.(a:t+„)^ A^A:j(T)fIri(a;t). 


0<(7<S 


(iii) Reaction fires and is not interrupted: 

Pimi^t — = 0}o<(T<T, p=i,...,Ai) = 1) 

X n Pi^Z,t = = 0}o<p<.. .=1,....A. Am = 1) 


0<(7<r 


X Piklt = lAi) 

(^1 - A/i,^(£Ct+tT)) A^A:i(T)fIri(a;(). 


(B2) 

(B3) 

(B4) 

(B5) 

(B6) 

(B7) 


(B8) 


(B9) 


0<cr<T 

pL—l,...^Ai 
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Following the same procedure as above and keeping track of the additional index /i, the generating functional 
is the same as when there is only one interruption reaction, only with the action 


= J dt 


p{t) ■ x{t) 


+ ^ I ^ dr e” Io - r)) 

/.=!,....A,-^0 ^ '' 

X J dTKi{T)nri{x{t — s))\ . (BIO) 
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